library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-6
library(lattice)
data(BCI)
?BCI
## UTM Coordinates (in metres)
BCI_xy = data.frame(x = rep(seq(625754, 626654, by=100), each=5), 
                    y = rep(seq(1011569,  1011969, by=100), len=50))
plot(BCI_xy)

Overall Data (just wanted to see)

sr = apply(BCI, 1, function(x) sum(x > 0))
hist(sr)

plot(BCI_xy, cex = sr/max(sr))

col_brks = hist(sr, plot=F)$breaks
col_indices = as.numeric(cut(sr, col_brks))
cols = rev(terrain.colors(length(col_brks)))
plot(BCI_xy, cex=2, pch=19, col=cols[col_indices])

# calculate Euclidean distance between richness and spatial coordinates
sr_dist = dist(sr)
xy_dist = dist(BCI_xy)
max_dist = max(xy_dist) / 2

# plot result
plot(xy_dist, sr_dist)
abline(lm(sr_dist ~ xy_dist), lwd=3, col='red')
lines(lowess(xy_dist, sr_dist), lwd=3, col='pink')
abline(v = max_dist, col='red', lwd=3, lty=2)

# compute correlation
obs_cor = cor(xy_dist, sr_dist)
obs_cor
## [1] 0.004339913
# carry out a permutation test for significance:
nperm = 1000
null_cor = obs_cor
for (i in 2:nperm) {
    # shuffle the rows of the spatial coordinates
    tmp_xy = BCI_xy[sample(nrow(BCI_xy)), ]
    # correlation between the shuffled spatial coordinates and sr_dist
    null_cor[i] = cor(dist(tmp_xy), sr_dist)
}
# compute the p-value
sum(null_cor >= obs_cor) / nperm 
## [1] 0.48
# carry out the same analysis
sr_mantel = mantel(xy_dist, sr_dist)
sr_mantel
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = xy_dist, ydis = sr_dist) 
## 
## Mantel statistic r: 0.00434 
##       Significance: 0.434 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0678 0.0883 0.1087 0.1290 
## Permutation: free
## Number of permutations: 999
  1. Common
alseis = BCI[,"Alseis.blackiana"]
alseis
##  [1] 25 26 18 23 16 14 18 14 16 14 14 19  8 17 15 25 31  7 13 10 12 22  5 14 20
## [26]  7 17 16 15 36 11 21 24 42 93  8 19 25 38 65 13 13  8 13 10 29 17 12  6  9
plot(BCI_xy, cex = alseis/max(alseis))

##ggplot solution
library(ggplot2)
ggplot(BCI_xy, aes(x=x, y=y, size=alseis)) + geom_point() + scale_size_continuous(range = c(1,5))

col_brks = hist(alseis, plot=F)$breaks
col_indices = as.numeric(cut(alseis, col_brks))
cols = rev(terrain.colors(length(col_brks)))
plot(BCI_xy, cex=2, pch=19, col=cols[col_indices])

# calculate Euclidean distance between abundance and spatial coordinates
ab_dist = dist(alseis)
xy_dist = dist(BCI_xy)
max_dist = max(xy_dist) / 2
plot(xy_dist, ab_dist)
abline(lm(ab_dist ~ xy_dist), lwd=3, col='red')
lines(lowess(xy_dist, ab_dist), lwd=3, col='pink')
abline(v = max_dist, col='red', lwd=3, lty=2)

# compute correlation
obs_cor = cor(xy_dist, ab_dist)
obs_cor
## [1] -0.02348007
# carry out a permutation test for significance:
nperm = 1000
null_cor = obs_cor
for (i in 2:nperm) {
    # shuffle the rows of the spatial coordinates
    tmp_xy = BCI_xy[sample(nrow(BCI_xy)), ]
    # correlation between the shuffled spatial coordinates and ab_dist
    null_cor[i] = cor(dist(tmp_xy), ab_dist)
}
# compute the p-value
sum(null_cor >= obs_cor) / nperm 
## [1] 0.593
ab_mantel = mantel(xy_dist, ab_dist)
ab_mantel
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = xy_dist, ydis = ab_dist) 
## 
## Mantel statistic r: -0.02348 
##       Significance: 0.603 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0818 0.1057 0.1210 0.1390 
## Permutation: free
## Number of permutations: 999
##examine community
## compute bray curtis distance for the community matrix
comm_dist = vegdist(BCI)
plot(xy_dist, comm_dist)
abline(lm(comm_dist ~ xy_dist), lwd=3, col='red')
lines(lowess(xy_dist, comm_dist), lwd=3, col='pink')
lines(lowess(xy_dist, comm_dist, f=0.1), lwd=3, col='blue')

abline(v = max_dist, col='red', lwd=3, lty=2)

comm_mantel = mantel(xy_dist, comm_dist)
comm_mantel
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = xy_dist, ydis = comm_dist) 
## 
## Mantel statistic r: 0.4078 
##       Significance: 0.001 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0642 0.0878 0.1019 0.1130 
## Permutation: free
## Number of permutations: 999
ab_corlog = mantel.correlog(ab_dist, xy_dist)
comm_corlog = mantel.correlog(comm_dist, xy_dist)
ab_corlog
## 
## Mantel Correlogram Analysis
## 
## Call:
##  
## mantel.correlog(D.eco = ab_dist, D.geo = xy_dist) 
## 
##         class.index     n.dist Mantel.cor Pr(Mantel) Pr(corrected)  
## D.cl.1   136.870241 144.000000   0.026911      0.105         0.105  
## D.cl.2   210.610723 376.000000  -0.015646      0.282         0.282  
## D.cl.3   284.351204 390.000000  -0.046374      0.028         0.084 .
## D.cl.4   358.091686 148.000000  -0.015208      0.230         0.460  
## D.cl.5   431.832168 372.000000  -0.044131      0.074         0.296  
## D.cl.6   505.572649 266.000000  -0.011319      0.179         0.537  
## D.cl.7   579.313131 168.000000         NA         NA            NA  
## D.cl.8   653.053613 100.000000         NA         NA            NA  
## D.cl.9   726.794094 154.000000         NA         NA            NA  
## D.cl.10  800.534576  88.000000         NA         NA            NA  
## D.cl.11  874.275058  50.000000         NA         NA            NA  
## D.cl.12  948.015539  24.000000         NA         NA            NA  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
comm_corlog
## 
## Mantel Correlogram Analysis
## 
## Call:
##  
## mantel.correlog(D.eco = comm_dist, D.geo = xy_dist) 
## 
##         class.index     n.dist Mantel.cor Pr(Mantel) Pr(corrected)    
## D.cl.1   136.870241 144.000000   0.181547      0.001         0.001 ***
## D.cl.2   210.610723 376.000000   0.133245      0.001         0.002 ** 
## D.cl.3   284.351204 390.000000  -0.017807      0.253         0.253    
## D.cl.4   358.091686 148.000000  -0.029965      0.091         0.182    
## D.cl.5   431.832168 372.000000  -0.055778      0.039         0.117    
## D.cl.6   505.572649 266.000000  -0.087321      0.001         0.006 ** 
## D.cl.7   579.313131 168.000000         NA         NA            NA    
## D.cl.8   653.053613 100.000000         NA         NA            NA    
## D.cl.9   726.794094 154.000000         NA         NA            NA    
## D.cl.10  800.534576  88.000000         NA         NA            NA    
## D.cl.11  874.275058  50.000000         NA         NA            NA    
## D.cl.12  948.015539  24.000000         NA         NA            NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##no significance in the abundance of Alseis but there is significant in the community structure (non linear)
par(mfrow=c(1,2))
plot(ab_corlog)
mtext(side=3, 'Species Abundance')
abline(v = max_dist, col='red', lwd=3, lty=2)
plot(comm_corlog)
mtext(side=3, 'Community Composition')
abline(v = max_dist, col='red', lwd=3, lty=2)

Rare

andira = BCI[,"Andira.inermis"]
andira
##  [1] 0 0 0 0 1 1 0 0 1 0 1 0 0 1 0 0 0 0 2 0 0 1 1 2 0 2 0 1 1 2 0 0 2 1 0 0 0 0
## [39] 0 0 1 1 2 0 0 0 0 1 2 1
plot(BCI_xy, cex = andira/max(andira))

##ggplot
ggplot(BCI_xy, aes(x=x, y=y, size=andira)) + geom_point() + scale_size_continuous(range = c(1,5))

col_brks = hist(andira, plot=F)$breaks
col_indices = as.numeric(cut(andira, col_brks))
cols = rev(terrain.colors(length(col_brks)))
plot(BCI_xy, cex=2, pch=19, col=cols[col_indices])

andira_dist = dist(andira)
xy_dist = dist(BCI_xy)
max_dist = max(xy_dist) / 2
plot(xy_dist, andira_dist)
abline(lm(andira_dist ~ xy_dist), lwd=3, col='red')
lines(lowess(xy_dist, andira_dist), lwd=3, col='pink')
abline(v = max_dist, col='red', lwd=3, lty=2)

andira_mantel = mantel(xy_dist, andira_dist)
andira_mantel
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = xy_dist, ydis = andira_dist) 
## 
## Mantel statistic r: -0.05389 
##       Significance: 0.889 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0640 0.0845 0.0982 0.1199 
## Permutation: free
## Number of permutations: 999
andira_corlog = mantel.correlog(andira_dist, xy_dist)
andira_corlog
## 
## Mantel Correlogram Analysis
## 
## Call:
##  
## mantel.correlog(D.eco = andira_dist, D.geo = xy_dist) 
## 
##         class.index      n.dist  Mantel.cor Pr(Mantel) Pr(corrected)
## D.cl.1  136.8702408 144.0000000  -0.0216764      0.233         0.233
## D.cl.2  210.6107225 376.0000000  -0.0324079      0.138         0.276
## D.cl.3  284.3512042 390.0000000  -0.0251373      0.201         0.414
## D.cl.4  358.0916859 148.0000000  -0.0143570      0.314         0.603
## D.cl.5  431.8321676 372.0000000  -0.0026076      0.486         0.804
## D.cl.6  505.5726492 266.0000000   0.0029042      0.499         1.000
## D.cl.7  579.3131309 168.0000000          NA         NA            NA
## D.cl.8  653.0536126 100.0000000          NA         NA            NA
## D.cl.9  726.7940943 154.0000000          NA         NA            NA
## D.cl.10 800.5345760  88.0000000          NA         NA            NA
## D.cl.11 874.2750577  50.0000000          NA         NA            NA
## D.cl.12 948.0155393  24.0000000          NA         NA            NA
par(mfrow=c(1,2))
plot(andira_corlog)
mtext(side=3, 'Species Abundance')
abline(v = max_dist, col='red', lwd=3, lty=2)

There is no evidence of spatial dependance in both the common (Alseis) and rare (Andira) species.

2.one species

sp_ids = c("Cordia.lasiocalyx", "Hirtella.triandra",
           "Picramnia.latifolia", "Quassia.amara",
           "Tabernaemontana.arborea", "Trattinnickia.aspera", 
           "Xylopia.macrantha")
BCI_ids = cbind("Drypetes.standleyi", "Cordia.lasiocalyx", "Hirtella.triandra",
           "Picramnia.latifolia", "Quassia.amara",
           "Tabernaemontana.arborea", "Trattinnickia.aspera", 
           "Xylopia.macrantha")
sp_a <- apply(BCI_ids, 1, function(x) sum(x > 0))
sp_data = data.frame(sp_a, BCI, BCI_xy)
BCI_sub = subset(sp_data, select = BCI_ids)
library(nlme)
sp_lm = gls(Drypetes.standleyi ~ Hirtella.triandra, data=BCI_sub)
var_sp = Variogram(sp_lm, form= ~ BCI_xy$x + BCI_xy$y)
plot(var_sp)

res = residuals(sp_lm)
plot(dist(sp_data[, c('x', 'y')]), dist(res))
lines(lowess(dist(sp_data[, c('x', 'y')]), dist(res)), col='red', lwd=2)
abline(v = max_dist, col='red', lwd=3, lty=2)

x = BCI_xy$x
y = BCI_xy$y
sp_exp = update(sp_lm, corr=corExp(form=~x + y))
plot(Variogram(sp_exp, maxDist = max_dist))

plot(Variogram(sp_exp, resType='normalized', maxDist = max_dist))

sp_exp_nug = update(sp_exp, corr=corExp(form = ~x +y , nugget=T))
plot(Variogram(sp_exp_nug, maxDist = max_dist))

sp_rat_nug = update(sp_lm, corr=corRatio(form = ~x + y, nugget=T))
plot(Variogram(sp_rat_nug, maxDist = max_dist))

plot(Variogram(sp_rat_nug, resType='n', maxDist = max_dist))

col_brks = hist(residuals(sp_rat_nug), plot=F)$breaks
col_indices = as.numeric(cut(residuals(sp_rat_nug), col_brks))
cols = rev(terrain.colors(length(col_brks)))
plot(BCI_xy, cex=2, pch=19, col=cols[col_indices])

anova(sp_lm, sp_exp, sp_exp_nug, sp_rat_nug, test=F)
##            Model df      AIC      BIC    logLik
## sp_lm          1  3 346.1696 351.7832 -170.0848
## sp_exp         2  4 312.5365 320.0213 -152.2682
## sp_exp_nug     3  5 313.2119 322.5679 -151.6059
## sp_rat_nug     4  5 310.6244 319.9804 -150.3122
summary(sp_lm)
## Generalized least squares fit by REML
##   Model: Drypetes.standleyi ~ Hirtella.triandra 
##   Data: BCI_sub 
##        AIC      BIC    logLik
##   346.1696 351.7832 -170.0848
## 
## Coefficients:
##                       Value Std.Error  t-value p-value
## (Intercept)       0.6267619 1.7550545 0.357118  0.7226
## Hirtella.triandra 0.3724844 0.1038112 3.588094  0.0008
## 
##  Correlation: 
##                   (Intr)
## Hirtella.triandra -0.806
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.5837021 -0.5405355 -0.2411211  0.1461580  4.0034096 
## 
## Residual standard error: 7.352136 
## Degrees of freedom: 50 total; 48 residual
summary(sp_exp)
## Generalized least squares fit by REML
##   Model: Drypetes.standleyi ~ Hirtella.triandra 
##   Data: BCI_sub 
##        AIC      BIC    logLik
##   312.5365 320.0213 -152.2682
## 
## Correlation Structure: Exponential spatial correlation
##  Formula: ~x + y 
##  Parameter estimate(s):
##        range 
## 1.399944e+13 
## 
## Coefficients:
##                       Value Std.Error      t-value p-value
## (Intercept)       14.473746 1746021.3  0.000008290   1.000
## Hirtella.triandra -0.000646       0.1 -0.006275689   0.995
## 
##  Correlation: 
##                   (Intr)
## Hirtella.triandra 0     
## 
## Standardized residuals:
##           Min            Q1           Med            Q3           Max 
## -8.288079e-06 -8.284289e-06 -7.138550e-06 -3.843500e-06  1.405581e-05 
## 
## Residual standard error: 1746021 
## Degrees of freedom: 50 total; 48 residual
summary(sp_exp_nug)
## Generalized least squares fit by REML
##   Model: Drypetes.standleyi ~ Hirtella.triandra 
##   Data: BCI_sub 
##        AIC      BIC    logLik
##   313.2119 322.5679 -151.6059
## 
## Correlation Structure: Exponential spatial correlation
##  Formula: ~x + y 
##  Parameter estimate(s):
##        range       nugget 
## 3.029648e+08 1.815151e-07 
## 
## Coefficients:
##                       Value Std.Error    t-value p-value
## (Intercept)       12.254798  6300.289 0.00194512  0.9985
## Hirtella.triandra  0.017615     0.099 0.17747775  0.8599
## 
##  Correlation: 
##                   (Intr)
## Hirtella.triandra 0     
## 
## Standardized residuals:
##           Min            Q1           Med            Q3           Max 
## -0.0019870533 -0.0019562986 -0.0016402511 -0.0007709315  0.0041779703 
## 
## Residual standard error: 6300.295 
## Degrees of freedom: 50 total; 48 residual
summary(sp_rat_nug)
## Generalized least squares fit by REML
##   Model: Drypetes.standleyi ~ Hirtella.triandra 
##   Data: BCI_sub 
##        AIC      BIC    logLik
##   310.6244 319.9804 -150.3122
## 
## Correlation Structure: Rational quadratic spatial correlation
##  Formula: ~x + y 
##  Parameter estimate(s):
##        range       nugget 
## 752.28335795   0.04166127 
## 
## Coefficients:
##                       Value Std.Error   t-value p-value
## (Intercept)       14.088134 14.874137 0.9471564  0.3483
## Hirtella.triandra  0.011441  0.087993 0.1300210  0.8971
## 
##  Correlation: 
##                   (Intr)
## Hirtella.triandra -0.167
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -0.7428201 -0.7362643 -0.6323781 -0.3420863  1.2834076 
## 
## Residual standard error: 19.19677 
## Degrees of freedom: 50 total; 48 residual

Rational model with the nugget fit better then without spatial error terms

all species

sp_all_lm = gls(Drypetes.standleyi ~ Cordia.lasiocalyx + Hirtella.triandra + Picramnia.latifolia + Quassia.amara + Tabernaemontana.arborea + Trattinnickia.aspera + Xylopia.macrantha , data=sp_data)
plot(Variogram(sp_all_lm, form= ~ BCI_xy$x + BCI_xy$y))

res = residuals(sp_all_lm)
plot(dist(sp_data[, c('x', 'y')]), dist(res))
lines(lowess(dist(sp_data[, c('x', 'y')]), dist(res)), col='red', lwd=2)
abline(v = max_dist, col='red', lwd=3, lty=2)

sp_all_exp = update(sp_all_lm, corr=corExp(form=~x + y))
plot(Variogram(sp_all_exp, maxDist = max_dist))

plot(Variogram(sp_all_exp, resType='normalized', maxDist = max_dist))

sp_all_exp_nug = update(sp_all_exp, corr=corExp(form=~x + y, nugget=T))
plot(Variogram(sp_all_exp_nug, maxDist = max_dist))

sp_all_rat_nug = update(sp_all_lm, corr=corRatio(form=~x + y, nugget=T))
plot(Variogram(sp_all_rat_nug, maxDist = max_dist))

plot(Variogram(sp_all_rat_nug, resType='n', maxDist = max_dist))

col_brks = hist(residuals(sp_all_rat_nug), plot=F)$breaks
col_indices = as.numeric(cut(residuals(sp_all_rat_nug), col_brks))
cols = rev(terrain.colors(length(col_brks)))
plot(BCI_xy, cex=2, pch=19, col=cols[col_indices])

col_brks = hist(residuals(sp_all_exp_nug), plot=F)$breaks
col_indices = as.numeric(cut(residuals(sp_all_exp_nug), col_brks))
cols = rev(terrain.colors(length(col_brks)))
plot(BCI_xy, cex=2, pch=19, col=cols[col_indices])

anova(sp_all_lm, sp_all_exp, sp_all_exp_nug, sp_all_rat_nug, test=F)
##                Model df      AIC      BIC    logLik
## sp_all_lm          1  9 307.1163 322.7554 -144.5582
## sp_all_exp         2 10 301.6062 318.9829 -140.8031
## sp_all_exp_nug     3 11 301.9592 321.0735 -139.9796
## sp_all_rat_nug     4 11 303.1486 322.2630 -140.5743
summary(sp_all_lm)
## Generalized least squares fit by REML
##   Model: Drypetes.standleyi ~ Cordia.lasiocalyx + Hirtella.triandra +      Picramnia.latifolia + Quassia.amara + Tabernaemontana.arborea +      Trattinnickia.aspera + Xylopia.macrantha 
##   Data: sp_data 
##        AIC      BIC    logLik
##   307.1163 322.7554 -144.5582
## 
## Coefficients:
##                             Value Std.Error   t-value p-value
## (Intercept)             -1.051752 2.1175346 -0.496687  0.6220
## Cordia.lasiocalyx        0.428920 0.2039316  2.103255  0.0415
## Hirtella.triandra        0.122279 0.0802638  1.523462  0.1351
## Picramnia.latifolia      0.662259 0.6358905  1.041468  0.3036
## Quassia.amara            4.085661 2.2842770  1.788602  0.0809
## Tabernaemontana.arborea -0.249725 0.1491192 -1.674667  0.1014
## Trattinnickia.aspera     1.349323 0.7147412  1.887848  0.0660
## Xylopia.macrantha        0.548832 0.1468772  3.736672  0.0006
## 
##  Correlation: 
##                         (Intr) Crd.ls Hrtll. Pcrmn. Qss.mr Tbrnm. Trttn.
## Cordia.lasiocalyx       -0.618                                          
## Hirtella.triandra       -0.212 -0.354                                   
## Picramnia.latifolia      0.025 -0.019 -0.381                            
## Quassia.amara            0.163 -0.378  0.307 -0.302                     
## Tabernaemontana.arborea -0.708  0.245  0.163 -0.113  0.148              
## Trattinnickia.aspera    -0.139  0.187 -0.311  0.308 -0.708 -0.144       
## Xylopia.macrantha       -0.140 -0.125  0.156 -0.463  0.314  0.279 -0.294
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.87708765 -0.42701500 -0.04032793  0.23615609  3.38768871 
## 
## Residual standard error: 4.539713 
## Degrees of freedom: 50 total; 42 residual
summary(sp_all_exp)
## Generalized least squares fit by REML
##   Model: Drypetes.standleyi ~ Cordia.lasiocalyx + Hirtella.triandra +      Picramnia.latifolia + Quassia.amara + Tabernaemontana.arborea +      Trattinnickia.aspera + Xylopia.macrantha 
##   Data: sp_data 
##        AIC      BIC    logLik
##   301.6062 318.9829 -140.8031
## 
## Correlation Structure: Exponential spatial correlation
##  Formula: ~x + y 
##  Parameter estimate(s):
##    range 
## 480.0567 
## 
## Coefficients:
##                             Value Std.Error  t-value p-value
## (Intercept)             2.3485197  6.154919 0.381568  0.7047
## Cordia.lasiocalyx       0.1208390  0.179811 0.672033  0.5052
## Hirtella.triandra       0.0191759  0.098501 0.194677  0.8466
## Picramnia.latifolia     0.2014516  0.509196 0.395627  0.6944
## Quassia.amara           1.2792289  1.847570 0.692385  0.4925
## Tabernaemontana.arborea 0.0674943  0.133782 0.504511  0.6165
## Trattinnickia.aspera    1.8115374  0.525147 3.449582  0.0013
## Xylopia.macrantha       0.3388574  0.156874 2.160064  0.0365
## 
##  Correlation: 
##                         (Intr) Crd.ls Hrtll. Pcrmn. Qss.mr Tbrnm. Trttn.
## Cordia.lasiocalyx       -0.226                                          
## Hirtella.triandra       -0.309 -0.022                                   
## Picramnia.latifolia      0.045 -0.066 -0.369                            
## Quassia.amara           -0.059 -0.304  0.321 -0.142                     
## Tabernaemontana.arborea -0.240 -0.016  0.288 -0.221  0.112              
## Trattinnickia.aspera    -0.069  0.168 -0.237  0.212 -0.633 -0.041       
## Xylopia.macrantha       -0.056 -0.137 -0.063  0.109  0.290  0.102 -0.186
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.0051632 -0.5235683 -0.3176178  0.2208753  2.3746027 
## 
## Residual standard error: 8.628464 
## Degrees of freedom: 50 total; 42 residual
summary(sp_all_exp_nug)
## Generalized least squares fit by REML
##   Model: Drypetes.standleyi ~ Cordia.lasiocalyx + Hirtella.triandra +      Picramnia.latifolia + Quassia.amara + Tabernaemontana.arborea +      Trattinnickia.aspera + Xylopia.macrantha 
##   Data: sp_data 
##        AIC      BIC    logLik
##   301.9592 321.0735 -139.9796
## 
## Correlation Structure: Exponential spatial correlation
##  Formula: ~x + y 
##  Parameter estimate(s):
##        range       nugget 
## 5.424635e+07 2.451077e-06 
## 
## Coefficients:
##                              Value Std.Error   t-value p-value
## (Intercept)              3.0501045 1785.0597  0.001709  0.9986
## Cordia.lasiocalyx        0.1426666    0.1895  0.752752  0.4558
## Hirtella.triandra       -0.0017714    0.0904 -0.019600  0.9845
## Picramnia.latifolia      0.2863335    0.5274  0.542880  0.5901
## Quassia.amara            1.3263643    1.9368  0.684818  0.4972
## Tabernaemontana.arborea  0.0407530    0.1395  0.292085  0.7717
## Trattinnickia.aspera     1.8170749    0.5730  3.171303  0.0028
## Xylopia.macrantha        0.4086700    0.1537  2.659324  0.0110
## 
##  Correlation: 
##                         (Intr) Crd.ls Hrtll. Pcrmn. Qss.mr Tbrnm. Trttn.
## Cordia.lasiocalyx       -0.001                                          
## Hirtella.triandra       -0.001 -0.098                                   
## Picramnia.latifolia      0.000  0.017 -0.360                            
## Quassia.amara            0.000 -0.292  0.344 -0.193                     
## Tabernaemontana.arborea -0.001 -0.020  0.160 -0.197  0.088              
## Trattinnickia.aspera     0.000  0.165 -0.276  0.255 -0.655 -0.036       
## Xylopia.macrantha        0.000 -0.066 -0.037 -0.048  0.306  0.140 -0.183
## 
## Standardized residuals:
##           Min            Q1           Med            Q3           Max 
## -0.0049325435 -0.0029088009 -0.0020671675  0.0005570041  0.0103697005 
## 
## Residual standard error: 1785.068 
## Degrees of freedom: 50 total; 42 residual
summary(sp_all_rat_nug)
## Generalized least squares fit by REML
##   Model: Drypetes.standleyi ~ Cordia.lasiocalyx + Hirtella.triandra +      Picramnia.latifolia + Quassia.amara + Tabernaemontana.arborea +      Trattinnickia.aspera + Xylopia.macrantha 
##   Data: sp_data 
##        AIC     BIC    logLik
##   303.1486 322.263 -140.5743
## 
## Correlation Structure: Rational quadratic spatial correlation
##  Formula: ~x + y 
##  Parameter estimate(s):
##       range      nugget 
## 402.2077831   0.2194023 
## 
## Coefficients:
##                             Value Std.Error   t-value p-value
## (Intercept)             2.0306920  5.171732 0.3926522  0.6966
## Cordia.lasiocalyx       0.1508099  0.194940 0.7736210  0.4435
## Hirtella.triandra       0.0076692  0.091987 0.0833720  0.9340
## Picramnia.latifolia     0.2509289  0.539635 0.4649976  0.6443
## Quassia.amara           1.5049423  1.960799 0.7675147  0.4471
## Tabernaemontana.arborea 0.0322219  0.142012 0.2268964  0.8216
## Trattinnickia.aspera    1.7698936  0.583930 3.0310015  0.0042
## Xylopia.macrantha       0.4058061  0.161181 2.5177087  0.0157
## 
##  Correlation: 
##                         (Intr) Crd.ls Hrtll. Pcrmn. Qss.mr Tbrnm. Trttn.
## Cordia.lasiocalyx       -0.273                                          
## Hirtella.triandra       -0.272 -0.122                                   
## Picramnia.latifolia      0.017  0.038 -0.387                            
## Quassia.amara           -0.039 -0.304  0.337 -0.213                     
## Tabernaemontana.arborea -0.242 -0.029  0.166 -0.201  0.106              
## Trattinnickia.aspera    -0.090  0.163 -0.272  0.271 -0.646 -0.036       
## Xylopia.macrantha       -0.095 -0.055 -0.073 -0.035  0.295  0.143 -0.164
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.0645964 -0.5625847 -0.3564599  0.2563989  2.6599508 
## 
## Residual standard error: 7.348772 
## Degrees of freedom: 50 total; 42 residual

Spatial error terms did not effect the model

Adding spatial error terms (rational nugget) for the model with one species improved the model. Fitting error to the model that included all species improved the model, but adding nuggets had no effect. Overall, coefficients for both models had an impact, but this was larger for the model with only one species. Spatial error terms can control for spatial dependence between neighbors and therefore improve our model.